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Abstract 

The accuracy of the numerical solution of a fractional differential 
equation depends on the differentiability class of the solution. The 
derivatives of the solutions of fractional differential equations often 
have a singularity at the initial point, which may result in a lower ac¬ 
curacy of the numerical solutions. We propose a method for improving 
the accuracy of the numerical solutions of the fractional relaxation and 
subdiffusion equations based on the fractional Taylor polynomials of 
the solution at the initial point. 
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1 Introduction 

Differential equations with fractional derivatives are used for modeling non- 
Markovian processes in various areas of science such as Physics, Biology 
and Economics [1-5]. The methods for analytical solution of ordinary and 
partial fractional differential equations can be applied to only special types 
of equations and initial conditions. In practical applications the solutions of 
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the fractional differential equations are determined by numerical methods. 
The development of efficient methods for numerical solution of fractional 
differential equations is an active research field. 

An important approximation for the Caputo fractional derivative is the 
LI approximation (l3lh It has been successfully used for numerical solution 
of ordinary and partial fractional differential equations. When the function 
y(x) is a smooth function, the accuracy of the LI approximation is O ( h 2 ~ a ). 
In previous work [l3| we determined the second order approximation ([4]) 
by modifying the first three coefficients of the LI approximation with the 
value of the Riemann zeta function at the point a — 1, and we computed the 
numerical solutions of the fractional relaxation and subdiffusion equations 
with sufficiently differentiable solutions. The solutions of fractional differen¬ 
tial equations often have a singularity at the initial point. In this case the 
numerical solutions may converge to the exact solution but their accuracy 
may be lower than the expected accuracy from the numerical analysis. In 
this paper we propose a method for improving the accuracy of the numerical 
solutions of the fractional relaxation equation 


y [a \x) + By(x) = 0, y( 0) = 1, 
where 0 < a < 1, and the fractional subdiffusion equation 
d a u(x,t ) d 2 u(x,t ) r . 

u(x , 0) = sin x, u( 0, t ) = 0, u( tt, t ) = 0. 


( 1 ) 

( 2 ) 


When the solutions of the relaxation and subdiffusion equations are smooth 
functions the accuracy of the numerical solutions using the LI approxima¬ 
tion are O (h 2 ~ a ) and O (t 2_q + h 2 ) and the numerical solutions using the 
modified LI approximation have accuracy O ( h 2 ) and O (t 2 + h 2 ). The solu¬ 
tions of the fractional relaxation and subdiffusion equations ([I]) and (j2J) have 
differentiable singularities at the initial point. The accuracy of the numerical 
solutions m an d (TTTT) of the fractional relaxation equation is O ( h a ) and the 
numerical solutions (I2U1) and (I2TT) of the fractional subdiffusion equation have 
accuracy O (r + h 2 ). Equations ([T]) and (j2J) are linear fractional differential 
equations. The form of the equations allows us to compute the Miller-Ross 
derivatives of the solutions at the initial point. In section 3 and section 
4 we use the fractional Taylor polynomials to improve the differentiability 
properties of the solutions and the accuracy of the numerical methods. 
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The Caputo derivative of order a, where 0 < a < 1 is defined as 


D%y(x) = y 


-- (Q) (x) = 


d a 

dx‘ 


;y(x) = 


r(i 


a 


m 

z-0' 


-df. 

V 3 


When the function y is defined on the interval (—oo, x], the lower limit of the 
integral in the definition of Caputo derivative is — oo. The value of the frac¬ 
tional derivative y <K ° i \x) depends on the values of the function on the whole 
interval of definition. One approach for discretizing the Caputo derivative 
is to divide the interval to subintervals of small length and approximate the 
derivative of the function y(x) using spline interpolation or a Lagrange poly¬ 


nomial 17]. Let h be a small positive number. Denote 


x n = nh , y n = y(x n ) = y{nh). 

The LI approximation for Caputo derivative is derived [131] by approximating 
the derivative of the function y{x) on the interval [xk-\,Xk\ with its value 
y'k-o .5 a t the midpoint of the interval. 


7/0 ~ 

an 


1 

h a T(2 


a) 


n— 1 

E (a) 

a k Vn-k , 

k =0 


(3) 


where a^ 0 = 1, = (n — l) 1 a — n 1 “ and 

4 a) = (k + l) 1 "" - 2k l ~ a + (k- l) 1_a , {k = 1, • ■ • , n - 1). 


When the function y(x) has a continuous second derivative, the accuracy of 
the LI approximation is O ( h 2 ~ a ) ([li|). The modified LI approximation for 
the Caputo derivative 


% 


(a) 


T(2 — a)h c 


X^4 a) Z/n-fc, 


(4) 


k =0 


has coefficients = cr^ for 2 < k < n, 

<5o a) = °‘o“ ) - co -1), 4 Q) = °i a) + 2 C(« - 1 ), 4 a) = 4 a) - CO - 1 ), 


where £(x) is the Riemann zeta function. The modified LI approximation 
has accuracy O (h 2 ). 
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By approximating the Caputo derivative at the point x n with Q and ([3]) 
we obtain the numerical solutions m and m of the fractional relaxation 
equation. In Table 1 and Table 2 we compute the maximum error and the 
order of the numerical solutions (TTOT) and (HI]) for the following equations 

y {0 ' 5) (x) + y(x) =x 2 + 2/(0) = 0, (5) 

oy 7T 

y (0 ' 5) (x ) + y(x) = x L25 + ^T 2 (0.25)x a75 , y( 0) = 0. (6) 

Equations (J5]) and © have solutions y(t) = x 2 and y(t) = x L25 . The solu¬ 
tion of equation (EJ) has an unbounded second derivative at t — 0. While 
the numerical solutions of equation © converge to the exact solution, their 
accuracy is smaller than O (h 2 ~ a ). 


Table 1: Maximum error and order of numerical solutions (TTOl) and (ITT]) for 


equation ([3]) 

on the interval [0,1]. 




h 

Error 

Order 

h 

Error 

Order 

0.05 

0.00281532 

1.45619 

0.05 

0.000695507 

1.90440 

0.025 

0.00101549 

1.47112 

0.025 

0.000182729 

1.92836 

0.0125 

0.00036392 

1.48049 

0.0125 

0.000047388 

1.94711 

0.00625 

0.00012986 

1.48661 

0.00625 

0.000012168 

1.96139 

0.003125 

0.00004621 

1.49072 

0.003125 

3.10 x 10~ 6 

1.97206 


Table 2: Maximum error and order of 
equation (J6J) on the interval [0,1]. 

numerical solutions (fTOl) and (TTTj) for 

h 

Error 

Order 

h 

Error 

Order 

0.05 

0.001825790 

1.1544 

0.05 

0.001825790 

1.15440 

0.025 

0.000806728 

1.17836 

0.025 

0.000806728 

1.17836 

0.0125 

0.000357787 

1.17298 

0.0125 

0.000351853 

1.19711 

0.00625 

0.000156707 

1.19103 

0.00625 

0.000151948 

1.21139 

0.003125 

0.000067873 

1.20717 

0.003125 

0.000065135 

1.22206 
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The fractional relaxation equation (jTO and the fractional subdiffusion 
equation ([2]) have solutions 

y{x) = E a (—Bx a ), u(x,t) = sinx-E a (— t a ). 


Both solutions have a differentiable singularity at the initial point of frac¬ 
tional differentiation. Jin, Lazarov and Zhou [16] discuss the numerical solu¬ 
tion of the fractional subdiffusion equation with accuracy in the time direc¬ 
tion O(r) and determine that the accuracy of the numerical solution (TTOl) of 


the relaxation equation (JTJ) is O ( h a ). Murillo and Yuste 22] use a finite dif¬ 


ference method with adaptive non-uniform timesteps for numerical solution 
of the fractional subdiffusion and diffusion-wave equations. In section 3 and 
section 4 we present a method for improving the accuracy of the numerical 
solutions of equations (H]) and (J2]) based on the fractional Taylor polynomials 
of the solution. 


2 Preliminaries 

The values of a differentiable function close to the point x = Xq are approxi¬ 
mated by its Taylor polynomials of degree m 

y (x) K 

n\ 

n= 0 

The definition of fractional Taylor polynomials at the point x = 0 is based on 
the composition of Caputo derivatives. The Millcr-Ross sequential fractional 
derivative for the Caputo derivative is defined as 

D ai y(x) = D%y{x), D ai+ol2 y(x) = D$D%y{x), 

and 

D ai+a2+ - +an y(x) = D^D* 2 ■ ■ ■ D^y(x). 

When y(x) is a differentiable function 

D°- 5+0 - 5 y(x) = D° c 5 y(x)D° c 5 y(x ) = D^y (x) = y'(x). 

The Caputo derivative D^y(x) = y'(x) and the Miller-Ross derivative D 0 5+0 ' 5 y(x ) 
of the function y(x) = 1 + x 0 ' 5 + x are different. 

D l c y(x ) = y'ix) = 1 + 0.5x -0 ' 5 , 
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D°™*y(*) = Dq 5 D Q 5 y(x) = D°S (r(1.5) + = 1- 

Two important functions in Fractional Calculus are the Mittag-Leffler and 
the Gamma function 

poo 

T(x) = / e~^ x ~ l di. 

Jo 

The Gamma function has values T(0) = T(l) = 1,T(0.5) = \fT and satisfies 

r(l + x) — aTLc), r(x)T(l — x) = —-—. 

sin ttx 

The one-parameter and two-parameter Mittag-Leffler functions are defined 
for a > 0 as 


EJx) = 


OO 


X n 


n r (an + 1) 

n =0 7 


> E a ^(x) — 'y ^ 


X 


T(an + /3) 

n —0 v 7 


When a — /3 — 1 the Mittag-Leffler functions are equal to the exponential 
function. The Mittag-Leffler functions appear in the solutions of fractional 
differential equations. The Laplace transform of the derivatives of the Mittag- 
Leffler functions satisfy 

Clx^E^ax-)} W = GLTM, 

(S a =F a) 

The next theorem is a generalization of the Mean-Value Theorem for differ¬ 
entiable functions. 


Theorem. (Generalized Mean-Value Theorem) 

y(x) = y( 0) + X y {a) (£ x ) (0 < ^ < x). 

1 (a + 1 ) 

In this paper we denote 

D na y(x) = D a+a+ - +a y(x ) = D a c D a c • • • D£y{x). 

When the function y(x) has nonzero Miller-Ross derivatives D na y( 0 ) the 
value of y(h ) is approximated by the fractional Taylor polynomials [10|. 


y(h) ~ 

n =0 


h an 

r (an + 1) 


D na y( 0), 
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where h is a small positive number. The Miller-Ross derivatives of the func¬ 
tion E a ( x a ) satisfy 

D na E a (: x a ) = E a 0 x a ), I7 nQ L; a (x 0 )!^ = (0) = 1. 

The integer order derivatives of the function E a ( x a ) —> ±oo as x —» 0. When 
h is a small number, the value of E a ( h a ) is approximated by its fractional 
Taylor polynomials. 

The fractional relaxation and subdiffusion equations, 

y {0 ' 5 \x) + y(x) = 0, 2/(0) = 1, (7) 

( d 05 u(x,t ) d 2 u(x,t) 

\ dt 05 = dx 2 ’ (8) 

i u(x, 0) = sin x , w(0, f) = u(7T, f) = 0, 

have solutions t/(x) = E 0 . 5 (—a; 0 ' 5 ) and u(x,t) = sin xE ( i 5 (—t 0 - 5 ). 



Figure 1: Graphs of the solutions of the fractional relaxation and subdiffusion 
equations (J7J) and (jSJ). 


In 13] we determined the hnite-difference approximations (1201) and (|2T]l 


for the fractional subdiffusion equation. The second derivative in the space 
direction is approximated by second order central difference. Numerical so¬ 
lution (T20l) uses the LI approximation for the Caputo derivative in the time 
direction and (l2Tj) uses the modified LI approximation. In Table 3 and Table 
4 we compute the maximal error and the order of numerical solutions (jTOj) 
and (ITT]) of the fractional relaxation equation (J7|) on the interval [0,1], and 
the numerical solutions ([20]) and d2T|) of the subdiffusion equation dHJ) on the 
rectangle [0,1] x [0,7r]. 
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Table 3: Maximum error and order of numerical solutions m and (1201) of 
the relaxation equation ([7|) on [0,1] (left) and the subdiffusion equation (JHJ) 
when r = 3h/7i, at time t = 1 (right). 


h 

Error 

Order 

r (h = nr/3) 

Error 

Order 

0.05 

0.0442319 

0.378959 

0.05 

0.00382782 

1.08802 

0.025 

0.0331978 

0.414001 

0.025 

0.00184202 

1.05524 

0.0125 

0.024484 

0.439247 

0.0125 

0.00089863 

1.03549 

0.00625 

0.0178342 

0.457192 

0.00625 

0.00044211 

1.02333 

0.003125 

0.0128769 

0.469859 

0.003125 

0.00021867 

1.01563 


Table 4: Maximum error and order of numerical solutions ca and (1211) of 
the relaxation equation ([7|) on [0,1] (left) and the subdiffusion equation (JH]) 
when r = 3h/7i, at time t = 1 (right). 


h 

Error 

Order 

r (h = nr/3) 

Error 

Order 

0.05 

0.0442319 

0.378959 

0.05 

0.00222875 

1.03303 

0.025 

0.0331978 

0.414001 

0.025 

0.00108596 

1.03727 

0.0125 

0.024484 

0.439247 

0.0125 

0.00053127 

1.03146 

0.00625 

0.0178342 

0.457192 

0.00625 

0.00026116 

1.02448 

0.003125 

0.0128769 

0.469859 

0.003125 

0.00012893 

1.01838 


The accuracy of numerical solution (TlOl) is O (h a ) ([16, Example 3.1]). 
The maximum error and order of numerical solutions m and CD for the 
fractional relaxation equation ([TJ) are given in Table 3 (left) and Table 4 (left). 
The two numerical solutions have the same maximal error. The maximum 
error is attained at the Erst approximation for y(h). Both numerical solutions 
use the same approximation V\ = w\ for y(h). Now we show that the error 
of the approximation V\ for y(h ) is O (h 0 5 ). 

00 (—h0-5\ n 

n=0 v ' 


,0.5 


7T 


+ 0{h). 


y(h ) = i 


2 























The approximation V\ = W\ for y[h ) is computed with v 0 = 1 and 


Vi - v 0 

T(1.5)/i° 


— + Vi — 0, V\ (l + r(1.5)/i°- J ) — 1 — 0, 


___= 1 _ h 0 - 5 + o (h) 

1 + T(1.5 )h 0 - 5 2 + 1 J 


Vi = 


\y(h) - Vl | = ( -y- - ) h °' 5 + 0 ( h ) ■ 


Then 


The accuracy of the approximation v\ = w\ for y(h) is O ( h °' 5 ). 


3 Numerical Solutions of the Fractional Relaxation Equation 


The fractional relaxation equation is an important ordinary fractional dif¬ 
ferential equation. The numerical and analytical solutions of the relaxation 
equation are discussed in [3,7,9,11-14]. The exact solution of the equation 

y ( °\x ) + By(x) = F(x), y( 0) = 1, (9) 


is determined with the Laplace transform method 

y(x) = E a (-Bx a ) + f X r-Xa (-BC) F(X - £)<%. 

Jo 

The numerical solutions {v n } and {w n } of the fractional relaxation equation 
(191) for approximations (J3J) and (j4]) are computed explicitly with v 0 = i(0). 


v n = 


+ Bh a T (2 - a) 
and w 0 — 1, wi — Vi, 

1 

8i a) + Bh a T {2 - a) 


T(2-a)h a F n -J2^^n-k\, ( 10 ) 


k =i 


w n = 


T(2 -a)h a F n -J2^^n-k). (11) 


k =i 


When F(x) = 0, the fractional relaxation equation 

y ( °\x ) + By(x) = 0, r/(0) = 1, 
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has the solution y(x) = E a (—Bx a ). In section 2 we computed the numerical 
solutions m and (fTTj) of equation (fill for a — 0.5. Both numerical solu¬ 
tions have maximum error at the first approximation v\ = w\ and accuracy 
O(h 0 5 ). The numerical experiments in this section confirm that the accu¬ 
racy of (flQj) and (TlTT) for equation (jTJ) is O (h a ) when a = 0.3 and a = 0.7 
and all values of a between 0 and 1. In this section we improve the accuracy 
of (fTUj) and (fTTjl using the fractional Taylor polynomials of the solution at 
the point x = 0. The improved numerical solutions have accuracy O (H 2 ~ a ) 
and O ( h 2 ). In Fig. 2 we compare (TTOl) with the improved numerical solution 
(TT4I) and the exact solution of equation (JTJ) for a = 0.3. 

The fractional relaxation equation is a linear ordinary fractional differen¬ 
tial equation. We can determine the Miller-Ross derivatives of the solution at 
x = 0 using fractional differentiation. The Miller-Ross fractional derivative 
of order a is equal to the Caputo derivative. Write the relaxation equation 
in the form 

D a y(x) + By(x) = 0, y(0) = 1. (12) 

When x = 0 we obtain D a y( 0) = —By{0) = —B. By applying fractional 
differentiation of order a 


D 2a y(x ) + BD a y(x ) = 0, 

D 3a y(x ) + BD 2a y(x ) = 0. 

From the above equations we determine the values of D 2a y{0 ) and D 3a y( 0). 

D 2a y( 0) = -BD a y{ 0) = B 2 , D 3a y( 0) = -BD 2a y{ 0) = —B 3 . 

By differentiating n — 1 times equation (fl2|) we obtain 

D na y( 0) = -BD {n - l)a y( 0). 


By induction we can show that D na y( 0) = (— B) n . When h is a small number, 
the fractional Taylor polynomials approximate the value of y(h) 


y(h) w 

n =0 


(—Bh a ) n 
r (an + 1) 


Now we use the fractional Taylor polynomials for improving the differen¬ 
tiability properties of the solution of the relaxation equation. Let m be a 
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positive integer such that ma > 2. Substitute 


z(x) = y(x) - 


y, ( -Bx a ) n 
^T(an+1)' 


The Caputo derivative of the function z(x) is 


m— 1 

z^\x) = y ( - a \x) + B s ^ 

n =0 


(■ -Bx a ) n 
T (cm + 1) 


The function z{x) is a twice continuously differentiable function and satisfies 
the relaxation equation 

/_ r>\m+l~*am 

z {a \x) + Bz(x) = p( am + t ) ■ *(0)=0. (13) 

The accuracy of numerical solutions (HU and (ITTh for equation (Tl3h is O ( h 2 ~ a ) 
and O ( h 2 ). Let {zk} be a numerical solution of ([T3]) . The improved numer¬ 
ical solution {y k } of equation (P) is determined from {z k } with 

(M) 

In Table 5, Table 6 and Table 7 we compute the numerical solutions of the 
relaxation equations (CD) and (1T3l) for a = 0.3 and a = 0.7. The number m is 
chosen such that ma > 2. 

• Let a = 0.3, B — 1 and m = 7. The fractional relaxation equation 

y {0 ' 3) {x) + y(x) = 0, y(0) = 1 (15) 

has analytical solution E 0 . 3 (—x 0 ' 3 ). Substitute 


z(x) =y( x ) + ^(-1) 


n+1 


X 


0.3n 


n=0 


T(0.3n + 1) ’ 


The function z(x) G C 2 [ 0,1] and is a solution of the equation 


£ (0 ' 3) (:r) +z(x) = 


x 


2.1 


r(3.i) 


, z(0) = 0. 


(16) 
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Figure 2: Graph of the analytical solution of equation dT5l) on the interval 
[0,1] and numerical solutions (HUjl -red and the improved numerical solution 
(U4]) -black for (fTOjh when h = 0.1. 



Table 5: Maximum error and order of numerical solutions (HUD and (TITT) for 
equation (TT5T) (left) and equation (H7D (right). 


h 

Error 

Order 

0.05 

0.0496688 

0.147528 

0.025 

0.0441259 

0.170715 

0.0125 

0.0386498 

0.191164 

0.00625 

0.0334386 

0.208945 

0.003125 

0.0286254 

0.224219 


h 

Error 

Order 

0.05 

0.0840166 

0.4486824 

0.025 

0.0567272 

0.566632 

0.0125 

0.0365217 

0.635288 

0.00625 

0.0229408 

0.670839 

0.003125 

0.0142437 

0.687593 


The accuracy of of numerical solutions (1T0|) and (ITT]) for equation (H5D is 
0(h) - Table 5 (left). The maximum error and order of numerical solutions 
(gSD and (fill) for equation (UUP are given in Table 6. In Figure 2 we compare 
the analytical solution of (113]) with numerical solution (HUD and its improved 
numerical solution (HU). 

• Let a = 0.7, m — 3 and B = 4. 

y^ 7 \x) + Ay(x) = 0, 2/(0) = 1. (17) 


The maximum error and order of numerical solutions (HUD and (HID for equa¬ 
tion (H7D are given in Table 5 (right). Substitute 


z(x ) = y(x) — 1 + 


4a; 1 


0.7 


16a; 1 ' 4 64a; 2 ' 1 

+ 


r(i.7) r(2.4) r(3.i)‘ 
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z(0) = 0. 


(18) 


The function z(x) is a solution of the equation 

r n7 ), , , s 256a; 2 - 1 

V ; y ' T(3.1) ’ 

The maximum error and order of numerical solutions (fTOl) and (HIT) for equa¬ 
tion ([HD are given in Table 7. 


Table 6: Maximum error and order of numerical solutions (fTOl) and (TTTh for 
equation (TT5h . 


h 

Error 

Order 

h 

Error 

Order 

0.05 

0.00022816 

1.61147 

0.05 

0.0000506703 

1.88353 

0.025 

0.00007356 

1.63301 

0.025 

0.0000132354 

1.93674 

0.0125 

0.00002347 

1.64829 

0.0125 

3.39 x 10~ 6 

1.96378 

0.00625 

7.43 x 10" 6 

1.65956 

0.00625 

8.61 x 10~ 7 

1.97867 

0.003125 

2.34 x 10" 6 

1.66807 

0.003125 

2.17 x 10~ 7 

1.98724 


Table 7; Maximum error and order of numerical solutions (fTUl) and (flip for 
equation (ITTjh 


h 

Error 

Order 

h 

Error 

Order 

0.05 

0.0825401 

1.27437 

0.05 

0.00682601 

2.21359 

0.025 

0.0338612 

1.28546 

0.025 

0.00136314 

2.32411 

0.0125 

0.0138328 

1.29154 

0.0125 

0.00024704 

2.46412 

0.00625 

0.0056374 

1.29499 

0.00625 

0.00004329 

2.51272 

0.003125 

0.0022943 

1.29701 

0.003125 

7.24 x 10~ 6 

2.57940 


4 Numerical Solutions of the Fractional Subdiffusion Equation 

The time-fractional subdiffusion equation is a parabolic partial fractional 
differential equation obtained from the heat equation by replacing the time 
derivative with a fractional derivative of order a, where 0 < a < 1, 

d a u(x, t) d 2 u(x,t ) , 1fVv 

8t a dx 2 ' ^ 


13 

























The analytical solution of the fractional subdiffusion equation with initial 
and boundary conditions 

u(x, 0) = Uo(x), u(0,t) = u(7r,t) = 0, 

is determined using the separation of variables method 

OO 

u(x,t) = c n E a (—n 2 t a ) sin(nx), 

n =1 


where 


c n = ~ / u 0 (£)sm(n£)d£,. 


7T 


Each term E a (—n 2 t a ) sin nx is a solution of (TT9|) and the coefficient c n is the 
coefficient of the Fourier sine series of the function u 0 (x). When u 0 (x) = sinx, 
the subdiffusion equation ([2]) has analytical solution u(x,t ) = sin xE a (—t a ). 
In section 3 we used the fractional Taylor polynomials of the solution of the 
relaxation equation to improve the accuracy of numerical solutions (fIU|) and 
([IT]) . In this section we use the fractional Taylor polynomials of the solution 
of the subdiffusion equation in the time direction for improving the accuracy 
of difference approximations (I20p and (121]) . 

Let h — tt/N,t — 1/M, where M and N are positive integers, and Q be 
a grid on the rectangle [0,7r] x [0,1] 


Q = {( nh , mr)||l < n < N, 1 < m < M} . 


In 13] we determined the finite difference approximations 
the the fractional subdiffusion equation 


and (l2Tf for 


d a u(x,t) d 2 u(x,t) 


+ F(x,t), 


dt a dx 2 

w(x,0) = uq(x), u(0,t) = u^n^t) = 0. 

Denote by u™ and F™ the values of the functions u(x,t) and F(x,t ) on Q. 
u™ = u(nh, mr), F™ = F(nh, mr). 


Let 


r] = r(2-a) —, 
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K be the tridiagonal matrix of dimension N — 1 with values 1 + 2 77 on the 
main diagonal, and —77 on the diagonals above and below the main diagonal 


K 5 


/l + 2r] -77 

—77 1 + 2rj 

0 —77 

0 0 

\ 0 0 


0 0 0 \ 

—77 0 0 

1 + 277 —77 0 

—77 1 + 277 —77 

0 —77 1 + 277 ) 


and V° be the vector of dimension IV—1 determined from the initial condition 


= Mn/OK • 

The difference approximation for the subdiffusion equation {V m } using the 
LI approximation (J3J) is computed with the linear system 

KV m = R u (20) 


where R\ is the vector 


Ri 


m 

-E I i" ) + + r " [ '(2 -“) 7 T 


The numerical solution {hh m } of the subdiffusion equation using the modified 
LI approximation (J3j) is computed with the linear system 


(K - C(a - l)/)l+ m = R 2 , 


( 21 ) 


for 2 < 777 < N — 1, where R 2 is the vector 


r 2 = 


k =1 


'w: 


+ t“T(2 — ot)F* 


-l N—l 


-I n= 1 


and W° = V°, W 1 = V" 1 . 

When the solution of the subdiffusion equation is a sufficiently differen¬ 
tiable function, the difference approximations (1201) and (12T| have accuracy 
O (r 2 ~ a + h 2 ) and O (r 2 + h 2 ). The solution of subdiffusion equation ([2]) has 
an unbounded first derivative at t — 0. The numerical experiments in section 
2 with a = 0.5 and other values of a G (0,1) indicate that the accuracy of 
numerical solutions (12Uj) and (Ell) is O (r + h 2 ). 
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Let D^u(x,t) be the Millcr-Ross derivative of order f3 of the function 
u(x,t) in the time direction. 


Dfu{x,t ) = 


d 2 u(x, t) 
dx 2 


Then 


D 2a u(x,t) = D^D^u(x,t) 


Df a u(x,t) = D?D 2a u(x,t) 
In this way we obtain 


n . a d 'Mx,t) 
t dx 2 

^ a d 4 u(x,t) 
t 8x 4 


d 2 

dx 2 


DfU(x, t) 


d 4 u(x, t) 
dx 4 


d 4 

dx 4 


D?u(x,t) 


d d u(x, t) 
dx 6 


Set t — 0 


D™u{x,t) 


d 2n u(x , t) 
dx 2n 


77rv . . d 2n u(x, 0 ) d 2n sin x . . 

= a ; 2 „ = = (—i) smz. 

From the fractional Taylor polynomials approximation 


u\x 


, h ) « sin x l) r 


h n 


T(na + 1) ’ 

n —0 v 7 

when h > 0 is a sufficiently small number. Substitute 

t na 


v(x,t) =u(x,t ) -sina: V(-l) n —-— 

“ 1 yna + 1) 

n =0 


We have that 


m—1 


AMM) = Dfu(x,t) +sinx y'(-l)”+ 1 ——- 

“ f yna + 1) 

n=0 

n2 n2 m ma 

u(x,f) = jj-^u(x, t) +sinx^(-l) n+1 


dx 2 


n =0 


T(na + 1) 
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The function v(x,t) is a solution of the fractional subdiffusion equation 


( d a v(x,t) _ d 2 v(x,t) 1 t ma 

l dt a dx 2 1 T(ma + 1 )’ 

[ v(x, 0 ) = u( 0 , t) = v(n, t ) = 0 . 

When a = 0.3 and m — 7, the fractional subdiffusion equation 


d 03 u(x,t ) d 2 u(x,t) 

dt 03 dx 2 ’ 

u(x, 0 ) = sin x, u( 0 , t) = u(n, t) = 0 , 


( 22 ) 


has analytical solution u(x,t ) = sinxT’o.s (—t 0 ' 3 )- The difference approxima¬ 
tions (1201) and (| 2 TT) for equation (122|) have accuracy 0{r) in the time direction 
(Table 8 (left) and Table 9 (left)). Let 


‘ J~j0.3n 

v(x,t) = u(x,t) — sinx> ( — l) n — - r. 

v ’ ’ v ’ ’ ^ ’ T(0.3n +1) 

The function v(x,t) is a twice continuously differentiable function and is a 
solution of the homogeneous subdiffusion equation 

f d°' 3 v(x,t) d 2 v(x,t) . t 21 

\ dt °- 3 = dx 2 + smx T{3J) , (23) 

[ v(x, 0) = v(0, t ) = v(tt, t ) = 0. 

In Table 8 (right) and Table 9 (right) we compute the maximum error and 
numerical order of difference approximations ( 120 ]) and () 2 T]) for the fractional 
subdiffusion equation (l23l) on the rectangle 

(x,t) e [o, tt] x [o, l], 

with step size in the space direction h = ttt/3. Numerical experiments with 
other values of the step sizes h and r confirm that the numerical solutions 
( 120 ]) and ( 12 TT) converge to the analytical solution of the fractional subdiffusion 
equation (l23l) with the expected accuracy O (r 2_ “ + h 2 ) and O (r 2 + h 2 ). A 
question for future work is to extend the method presented in this paper to 
other ordinary and partial fractional differential equations. 
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Table 8: Maximum error and order of numerical solution (1201) of equations 
(|22|) (left) and (1231) (right), for h = 7tt/ 3, at time t — 1. 


T 

Error 

Order 

0.05 

0.00208324 

1.08653 

0.025 

0.00100782 

1.04759 

0.0125 

0.00049480 

1.02633 

0.00625 

0.00024489 

1.01468 

0.003125 

0.00012175 

1.00826 


r 

Error 

Order 

0.05 

0.000247016 

1.64590 

0.025 

0.000078268 

1.65812 

0.0125 

0.000024643 

1.66722 

0.00625 

7.72 x 10~ 1 2 3 4 5 6 

1.67411 

0.003125 

2.41 x 10~ 6 

1.67940 


Table 9: Maximum error and order of numerical solution ()2TT) of equations 
(122|) (left) and (123|) (right), for h = ttt/3, at time t — 1. 


T 

Error 

Order 

0.05 

0.001565250 

1.05230 

0.025 

0.000760933 

1.04055 

0.0125 

0.000372634 

1.03001 

0.00625 

0.000183395 

1.02281 

0.003125 

0.000090563 

1.01796 


T 

Error 

Order 

0.05 

0.000031877 

1.81175 

0.025 

8.54 x 10~ 6 

1.90107 

0.0125 

2.22 x 10' 6 

1.94430 

0.00625 

5.67 x 10~ 7 

1.96751 

0.003125 

1.44 x 10' 7 

1.98068 
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